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I. INTRODUCTION 

Ultracold gases of atoms manipulated by magnetic and 
laser fields proved to be a powerful tool to study quan- 
tum many-body problems [1, 2]. Paradigmatic models, 
such as the Hubbard model, which have been introduced 
to phenomenologically describe fundamental effects of 
solid-state physics, can now be realized experimentally 
with ultracold atoms [3]. The atom-based implementa- 
tion also offers the intriguing possibility of tuning the sys- 
tem parameters and, thereby, reaching so-far unexplored 
regimes of these models. Thus, beyond quantum simula- 
tion of various many-body problems [4], the high-degree 
of control of the dynamics in a reduced Hilbert-space can 
lead to novel applications in quantum measurement and 
also in quantum information processing. 

Experimental setups of cavity quantum electrodynam- 
ics (cavity QED) allow for the realization of quantum 
gases with long-range interaction [5-9] . Instead of short- 
range collisions, the dominant atom-atom interaction is 
mediated by the electromagnetic radiation field enclosed 
into a high-finesse resonator [10, 11]. The many-body 
ensemble couples globally to the resonator field, which 
situation, in certain cases, can be described by collective 
spin-boson models. In particular, there has been recently 
a breakthrough in the realization of the Dicke model [12]. 
It is a paradigmatic model of quantum optics describing 
the interaction of a single mode of the radiation field with 
a collection of two- level atoms (spins), 

H = u c a) a + uir S z + —= + a) S x , (1) 

v N 

where a is the bosonic annihilation operator of the field 
mode, and the collective spin S = ~^2u = i^ % \ the index 
i labeling the spin <x associated with the individual two- 
level atoms. The main interest in this model is related to 
the predicted quantum phase transition [13-15] between 
the normal phase (all atom is in the ground state, and 
the radiation field is vacuum) and a superradiant phase 
(the atoms and the mode have a coherent polarization 

(^Sj 7^ and field amplitude (a) ^ 0, respectively) at 
the coupling strength y = Ju)c wr, i.e., when it reaches 



the geometric mean of the eigenfrequencies ujq an d ojr of 
the field mode and the atoms, respectively. However, it 
is generally believed impossible to achieve such a strong 
dipole coupling between atoms and the radiation field in 
practice. 

A completely different approach, based on the mo- 
tional excitations of a Bose-Einstein condensate loaded 
into an optical cavity, led to the first observation of the 
phase transition in the Dicke model [16]. The Hamilto- 
nian describing the coherent coupling of a single mode 
of the resonator with two momentum eigenstates of a 
degenerate quantum gas driven by a laser perpendicular 
to the cavity axis (see Fig. 1) was shown to be formally 
equivalent to the Dicke model [16, 17]. Beyond the map- 
ping of the phase diagram [16], the spontaneous symme- 
try breaking [18] as well as the mode softening in the 
excitations spectrum at the critical point [19] have been 
demonstrated experimentally. This approach does not in- 
volve electronic states of the atoms other than the ground 
state. The characteristic frequency scale of the interac- 
tion is determined by the recoil frequency lor = k 2 /2m 
of the specific atom (k is the wavenumber of the driving 
laser, m is the atomic mass), which is typically in the kHz 
range, well below the internal, electronic energy scale. 
The coupling strength can be tuned continuously by ad- 
justing the power of the driving laser field. The Dicke 
model phase transition with ultracold atoms is, in fact, 
the zero temperature limit of the atomic self-organization 
in a cavity, which has been predicted in [20] and demon- 
strated in [21] for cold but not ultracold atoms. 

The cavity-based realization is an open system. The 
atoms are driven by a monochromatic laser source, scat- 
ter photons into the cavity which then leak out through 
the end mirrors. Therefore the Dicke-type Hamiltonian 
does not provide for a full description of the system. The 
critical behavior has been reinvestigated recently for the 
stationary state of the driven and damped system. Using 
a generalized Bogoliubov theory to treat the thermody- 
namic limit, it was shown that the non-equilibrium sys- 
tem also exhibits a dynamical quantum (T — 0) phase 
transition [22]. The critical point as well as the mean 
field solution are only slightly modified with respect to 
the phase transition in the ground state of (1). However, 
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the correlation functions describing the quantum fluctua- 
tions differ significantly in the two, equilibrium and non- 
equilibrium, cases. This difference can be captured, for 
example, by looking at the critical exponents of the di- 
verging photon number close to the critical point. There 
is also a drastic modification of the atom-field entangle- 
ment, since the singularity of the logarithmic negativity 
of the ground state [23] at the critical point is regular- 
ized, though some finite entanglement is still found in the 
steady state of the driven open system. 

In this paper we study the finite-size scaling in the 
phase transition of the open-system Dicke model. On 
the one hand, as the experiments involve obviously a fi- 
nite number of atoms, it is necessary to calculate the 
correlation functions for a finite system. On the other 
hand, theoretically, the finite-size scaling gives insight 
into the criticality and its classification. Various systems 
can show similar critical behavior, that leads to the con- 
cept of universality: the systems sharing the same criti- 
cal exponents correspond to the same universality class. 
For the ground state of the Dicke model the finite-size 
corrections and scaling exponents are known [24], and 
show a good agreement with the critical scaling behav- 
ior of infinitely coordinated spin systems described by 
the Lipkin-Meskhov-Glick model [25]. This correspon- 
dence has been investigated both analytically [26] and 
by numerical calculations [27]. It is well known that for 
infinitely coordinated systems the mean-field approach 
is exact in the thermodynamic limit and the large-size 
critical behavior is related to the upper critical dimen- 
sionality of the corresponding short range system [28]. 
Note that the short-range equivalent of the Dicke model 
is unknown. Here we will consider exponents for a non- 
equilibrium quantum critical system, for the open Dicke 
model, in order to make a first step towards the classifi- 
cation. With the help of a scaling hypothesis consistent 
with the numerical findings we identify the two critical 
exponents describing the universal behavior of the sys- 
tem. 

The mean-field theory we used previously is suitable 
to describe the thermodynamic limit [22]. In this paper 
we develop a Hartree-Fock-Bogoliubov (HFB) theory for 
the finite open system. The population of the excitation 
modes above the condensate is not due to a finite temper- 
ature but the atom-light interaction leads to a quantum 
depletion of the ground state. Close to the critical point, 
a macroscopic part of the atoms can be out of the conden- 
sate. If the total number of atoms is finite, the quantum 
fluctuations represented by such a macroscopic depletion 
are expected to yield a considerable back action on the 
mean field. This effect is taken into account by the HFB 
theory. 

The structure of this paper is as follows. In Sec. II we 
present the system of ultracold atoms coupled to a cavity 
field together with the physical approximations adopted 
in the description, and introduce the basic equations of 
motion. In Sec. Ill we derive the HFB approximation 
generalized to an open system. The theory involves a 




FIG. 1: (Color online) Self-organization phase transition of a 
BEC in a cavity. Below a threshold in the transverse driving 
field (left) the condensate is quasi-homogeneous, and there is 
no radiation field inside the cavity. Above threshold (right), a 
standing matter wave of period A appears that scatter photons 
into the cavity. 

self-consistency loop, therefore, a better understanding 
of the approach can be gained by the details of the it- 
erative algorithm of the solution which we summarize in 
a separate subsection. The results are presented then 
in Sec. IV. The phase transition is exhibited, however, 
we focus mostly on the finite-size scaling of the quantum 
criticality. We conclude and give an outlook in Sec. V. 
Some lengthy expressions are showed in Appendices. 

II. DESCRIPTION OF THE SYSTEM 

We study the quantum phase transition of a Bose- 
Einstein condensate which is placed into an optical res- 
onator and pumped by a laser from a direction perpen- 
dicular to the resonator axis [29-32] . We consider a single 
cavity mode which is quasi-resonant with the laser pump 
frequency w. The detuning between the laser and the 
atom, defined as = uj — u>a, is large enough so that 
the internal dynamics of the atoms can be adiabatically 
eliminated and only the center-of-mass motion takes part 
in the dynamics. For simplicity, the motion is considered 
in one dimension along the cavity axis direction denoted 
by x. The cavity length is L and the mode function in 
the relevant direction is f(x) = cos(kx). The atoms are 
represented by the second-quantized boson field opera- 
tor \& (aj) , whereas the resonator mode is described by 
the operator a. 

The Hamiltonian in units of h — 1 is 

+ Uq a) a cos 2 (kx) 
+ r} t cos(fcr) (a f e~ iut + ae itJt )^J $(z) dx . (2) 

The first term gives the energy of cavity photons of fre- 
quency ujc, the second one the motional energy of the 
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atom field. This latter is characterized by the so-called 
recoil frequency ujr = fc 2 /2m. The dispersive atom-light 
interaction yields two terms. Firstly, there is an atomic 

a 2 

refractive index term with coupling strength Uq = 
where g is the single-photon Rabi frequency of the cavity 
mode. Secondly, there appears an effective cavity driv- 
ing by means of Raman scattering from the transverse 
pump laser with amplitude r\ t — wnere is 

the laser Rabi frequency. In this term the time evolu- 
tion at the pump frequency u> is made explicit and we 
used the rotating wave approximation. We neglected the 
atom-atom s-wave collision term, because the optical in- 
teraction dominates in the experiments. 

The system is spatially periodic with respect to the 
cavity mode wavelength, therefore a discrete Fourier ba- 
sis is convenient to decompose the atom field, 
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(3) 



Since parity is conserved in the system, no transition 
from the even cosine functions to the odd sine functions 
is allowed. Assuming initially a homogeneous conden- 
sate, we can keep only the subspace of the cosine modes. 
Let us introduce the notation c = (c , ci , c 2 , ... Y for 
a column, and = ^cj , c\ , c\ , ■■■^j for a row vector 

to make the forthcoming expressions compact. For ex- 
ample, the total number of atoms reads N = c and 

N = ^-^y Obviously, in the numerical calculation, the 
Fourier expansion (3) is truncated at a cutoff index n ma _ x . 

The Hamiltonian (2) conserves the atom number, i.e., 
H and N commute, which enables us to introduce the 
grand canonical Hamiltonian given by K = H — fiN, 
where [i is the chemical potential. The grand canonical 
Hamiltonian can be built up from quadratic forms, and 
reads in a frame rotating at the laser frequency uj 



K = -A c a f a + ojr (c f M (0) c 

(ot+o) (c'MWc 



Vt 

+ ~U Q a< a (c f M (2) c) -//c f c, (4a) 

where Ap = cj — wp. The kernel matrices are given 
by: 
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(4d) 



Note that the matrix is a real, symmetric band 

matrix and the left and right bandwidths are given by 
the index j. This means that the laser driving couples 
only adjacent modes, while the light shift term couples 
second neighbors. 

The open system description, taking into account the 
irreversible loss of photons through the out-coupling mir- 
ror, relies on the Heisenberg-Langevin equations of mo- 
tion [33] 



i — a(t)= a[t),K - in a{t) + i£(t) (5a) 



(5b) 



where the cavity photon loss rate is 2k, and £(t) is a white 
noise. The noise operator has zero mean value, and it's 
only non- vanishing second-order correlation function is 



(l(f)l t (0) = 2**(*-0- 



(6) 




We use mean-field approach to describe the system. 
The boson operators of the photon and atomic modes 
are split into mean-field amplitudes and fluctuations, 

(7a) 
(7b) 

where (a) = and (c) = by definition, a is the co- 
herent cavity field and 7 is the condensate wave-function 
in the basis given by (3). This latter vector is normalized 
as 7^ -7 = 1, thereby N c is the number of condensed 
atoms. We assume that the value of N c is a constant of 
time [41]. Finally, the operators a and c account for the 
quantum fluctuations around the coherent mean value, 
e.g. for the atoms outside the condensate. 

The coupling constants in the interaction terms can be 
redefined as 



y = y/2N~ c TH 

u=-N c U , 



(8a) 
(8b) 
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The new parameters u and y have a constant value in 
the thermodynamic limit which is defined by N — > oo 
and L — > oo, while keeping the density N/L constant. 
This follows from the fact that g ~ l/\/L, which leads 
to i] t ~ 1/VL and U ~ 1/L. 

On substituting the mean-field decomposition into the 
equations of motion (5), a hierarchy of terms according 
to the different orders of the a and c operators can be 
established. As this intermediate step is needed only for 
the derivation, we display the lengthy expressions in Ap- 
pendix A, in Eqs. (Al) and (A2). These include the 
equations of motion both for the mean fields and for the 
fluctuation operators. To make the set of equations com- 
plete, one needs the relation 



N r = N- 



c' c 



(9) 



for the number of condensed atoms, which follows from 
the mean-field decomposition and from the normaliza- 
tion of 7. Without further approximations, the set of 
equations involve nonlinear, operator valued differential 
equations, which cannot be solved in practice. 



III. HARTREE-FOCK-BOGOLIUBOV THEORY 

In this section we present how we go beyond the 
Bogoliubov-type mean field analysis that we used pre- 
viously to describe the damped-driven system of a cav- 
ity mode coupled to motional excitations of a BEC [22]. 
The method is a kind of Hartree-Fock-Bogoliubov (HFB) 
approximation [34] generalized to the steady-state of an 
open system. It is needed although the system is at 
T = 0, because the values of the second order correlation 
functions are large near the phase transition point. The 
higher-order terms, neglected in the bare Bogoliubov- 
approximation, can have a considerable effect on the 
mean value equations. Moreover, the HFB approach 
gives atom number dependent corrections to the results 
obtained for the thermodynamic limit, which we will use 
for the finite-size scaling. 



A. The self-consistent mean field equations 

The starting point of the theory is the set of Eqs. (Al) 
and (A2). In the simple Bogoliubov- approximation, the 
terms of second and of third order in the fluctuation op- 
erators are simply neglected, so we obviously get linear 
equations for the fluctuation operators. At variance, in 
the HFB method, these higher order terms are approxi- 
mated by using lower order terms as shown in Appendix 
A. The second order terms are substituted by their ex- 
pectation values. In the third order terms, we substitute 
the product of two operators by its expectation value, 
and leave the third one as an operator. There are three 
possible ways of doing this, so we do it in all possible 
ways and each third order term is approximated by the 



sum of three linear terms. The equations of motion for 
the fluctuation operators become then linear in the HFB 
approximation. 

The HFB equation of motion for the cavity mean field 

is 

i^-a={n-in)a+\y f 7 t M (1) 7+^- (c^c)) 
+ ~ (^c t 5^M (2) 7 + 7 t M (2) (co>) , (10a) 
and for the condensate it reads 

1 



yM w (<a f c) + (ac)) 



a*M (2) (ac) +aM (2) (5 f c) 



(10b) 



Here we defined the renormalized cavity frequency 

il = -A c + u 7 f M (2) 7 + ± ^M (2) c) , (11) 

which includes the resonance shift due to the dispersive 
interaction with the atoms. We also introduce the real 
and symmetric matrix 

M = u R M (0) + -y {a* + a) M (1) 

+ ua*aM^ + ^ (^d)M {2 \ (12) 

which includes the eigenenergy and the coupling between 
the atomic modes via the cavity field. This coupling will 
be diagonalized in the next subsection. 

The equation of motion for the fluctuating term of the 
cavity field is 



i^a = (n-in)a + ua ^7 f M (2) c + c f M (2) 7 
X -y ( 7 t M (1) c + c t M (1) 7 



+ 



u 



c f a) M (2) c + c f M (2) (co>)+ i£(t) . (13a) 



The dynamics of the atomic excitation modes above the 
condensate is 



i^c = (M - (J.) c + -y (a) + a) M (1) 7 



it 



+ u [a* a + a^a) M (2) 7 

M (2) (& c) + h M (2) (ate)] . (13b) 



Equations (10) and (13) fully define the HFB theory. 
The terms proportional to 1 /N c are new with respect to 
the Bogoliubov-approach. They express the back-action 
of the fluctuations on the mean values. Therefore the 
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mean field and the correlation functions have to be de- 
termined self-consistently. We will look for the steady 
state solution of the coupled mean-field and fluctuation 
equations. This means that the mean fields a and 7 as 
well as the second order correlation functions built from 
a(t) and c(t) have a constant value in time. As we noted 
earlier, it would be enough to set the condensate number 
N c a constant of time at this point. The fluctuation op- 
erators, however, still strongly fluctuate in time because 
of the quantum noise £(t) driving terms. 



C. Linearized equations of fluctuations 

The equations of the fluctuation operators in the new 
basis are 



i — a = (n-in)a + ua( /3 f M (23 b + b f M (2) /3 
at V 



+ 1 -y (/3tM (1) b>b f M (1) /3 



£ ((h f ~aW% + b f M (2) (ba 



)+if(i), (20a) 



B. Decoupling the atomic modes 

In the first step, we resolve the direct coupling be- 
tween the atomic modes c as described in Eq. (13b) by 
the matrix M. Since M is real and symmetric, it can be 
diagonalized by an orthogonal transformation O, where 
O O t = O t = 1. Let us perform the transformation 
on the atomic modes 7 = O j3 and c = O b and simulta- 
neously 



M 



O t MO 



M (J) -> O t M U) O = M 



,~t0') 



(14a) 
(14b) 



where the matrix A is diagonal by the definition of the 
orthogonal matrix O. Thanks to the orthogonality of the 
transformation, the vector (3 is also normalized: 0* P = 1 
and the number of condensed atoms is given by 



N r = N - 



b f b 



(15) 



The steady-state mean-field equations in the new basis 
read 



= (fi 



u 



■)a+ l -y ^M (1) /3+i- 
((b t 5)M (2) /3 + /3tM (2) 



b f M (1) b 



b a 



= (A - fit) p + R 

where f2 in the transformed basis is 

u 

N~c 

and the new quantity R was introduced: 



= -A c + u P^M (2) P- 



b f M (2) b 



(16a) 
(16b) 

(17) 



R 



1 



^M (1) ((«t b 



+ u [a 



(Vm (2) (ah) 



ah 



aM (2) / 5 f b 



(18) 



The value of the chemical potential can be determined 
by multiplying both sides of Eq. (16b) by /r from the 
left, 

n = pi Ap + p^n, (19) 

and we made use of the fact that P is normalized to one. 



i|b=(A-/xI)b+^(at + a) M W /3 

(2) 



It 



+ u (a* a + a f a) M P 
fit M (2) (5b)+aM (2) (fitb 



(20b) 



Note that, in the new basis, the b operators 
are also decoupled from each other. Let us put 
the fluctuation operators into vector form, v(t) = 
/ ~t \ T 

(a(t), aS(t), h(t), h (t)\ , then the coupled linear 
equations (20) can be written in the compact form 



ij t Ht) = Fv(t) + iq(t) 



(21) 



where the F matrix contains the coefficients and q(t) = 

T 



£(£), &(t), 0, 0) is the noise. The noise correlations 
are given by 



(q n (t)q j {t'))=D nj S(t-t') 
where the D diffusion matrix reads 



D 



2k 




(22) 



(23) 







We have to determine the eigenvalues and the left and 
right eigenvectors of the matrix F: 



F r (j-) = uij r (j) 
F f 1 0) = uo* . 



(24a) 
(24b) 



These vectors form a bi-orthogonal basis, 1^ t — 
5k j . The vector v(i) can be decomposed in terms of 
quasi-normal modes 



U) 



The equation of pk(t) reads 
d 



dt 



Pk(t) 



iu k pk(t) + Qk{t) 



(25) 



(26) 



G 



where Qt(t) = l^^q(i) is the projected noise. The solu- 
tion is readily obtained 

Pk(t) = p k (0)e- iu " t + f dt'Q k {t')e- 1 ^^ . (27) 



After some straightforward calculations, we get for the 
second order correlation functions in the steady state 



n, j k,l 



(Ot 



1 



i(u>k + u>{) 



r (i) 

' fj, ' V 



(28) 

These are the second-order correlation functions needed 
for the HFB theory both for the mean field and for the 
fluctuation equations. Let us emphasize at this point 
that the source of non- vanishing correlations is the quan- 
tum noise associated with the dissipative dynamics of the 
system. In this model the only dissipative process is the 
photon field decay, and correspondingly, the correlation 
functions are proportional to n. 



8. Use Eq. (28) to get the new second-order correla- 
tion functions. 

9. This step completes the iteration algorithm, return 
to step 1. 

The iteration is performed until it converges, which is 
tested, in our routine, by the variation of the solution for 

a. At the end, we have the condensate wave function /3 
and the second-order correlation functions of the modes 

b. To obtain them in the original basis, 7 and the cor- 
relation functions of the modes c, one needs to perform 
the inverse of the orthogonal transformation O, which 
has anyway been determined during the iteration. 



IV. FINITE-SIZE SCALING 

We showed that in the thermodynamic limit the open 
system exhibits critical behavior in the vicinity of the 
critical coupling [22] 



D. Iteration algorithm for the self-consistent 
solution 



Vc 



(29) 



In order to solve the coupled equations of the HFB 
theory self-consistently, we adopt an iteration algorithm. 
Here we present it in some detail because the algorithm 
itself sheds more light on the structure of the cross- 
coupled mean-field and fluctuation equations. 

The fixed system parameters are N, ojr, Ac, f]ti Uq 
and k. Initially, the a and (3 mean-fields, and all the 
second-order correlation functions receive a random ini- 
tial value. The iteration algorithm is the following: 

1. Determine the total number of condensed atoms 
N c from Eq. (15) and calculate then the coupling 
constants y and u from Eq. (8). 

2. Determine the matrix M from Eq. (12). 

3. Decouple the atomic modes: diagonalize the matrix 
M to get the orthogonal transformation O accord- 
ing to Eq. (14a), and then perform the transforma- 

~ (J) 

tion of Eq. (14b) to obtain the matrices M 

4. Calculate R from (18) and evaluate the chemical 
potential fj, from Eq. (19). 

5. Update the value of (3 by using Eq. (16b); note that 
the diagonal matrix (A — /xl) is straightforwardly 
inverted. 

6. Update the value of the mean-field a from Eq. (16a) 
by using Q from Eq. (17). 

7. Having the new mean field values, the coefficients 
of the linearized differential equation (20) can be 
inserted into the matrix F of (21), and the eigen- 
values, the left- and right-eigenvectors are to be 
calculated [421. 



where 8c = —A c + \NU is the shifted frequency of the 
cavity mode. In the normal phase (y < y c ), the conden- 
sate wave function is homogeneous 7 = (1,0,0,...) and 
the mean cavity field is zero ao = 0, while in the super- 
radiant phase (y > y c ) a finite mean field builds up in 
the higher momentum modes of the condensate (7^ 7^ 0) , 
and coherent light field appears in the cavity (|ao| > 0). 
On top of the mean fields, there are fluctuations strongly 
increasing near the critical point. The correlation func- 
tions of the fluctuations are singular and diverge in the 
Bogoliubov theory according to scaling laws. This non- 
equilibrium phase transition is accompanied by critical 
slowing down, which is responsible for the scaling behav- 
ior. Quite remarkably the order parameter scales in a 
similar way as for the closed system. However, the sus- 
ceptibilities characterizing the fluctuations are governed 
by completely different scaling laws. In this section we 
focus our attention mainly on the effects of finite particle 
number, that is finite-size scaling. 



A. Regularization of the criticality 

In the HFB theory, we couple the correlation functions 
of the fluctuations back to the mean field equations. For 
finite atom number, this correction regularizes the critical 
divergence of the fluctuations, and account for nontrivial 
finite size effects and scaling. Fig. 2 presents the cavity 
photon number and shows the convergence towards the 
singularity of the thermodynamic limit when the atom 
number is gradually increased. The photon number is 
split, according to the mean-field separation (7), 



(at a) = Nc | a |2 + ( g t 2 ) 



(30) 
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FIG. 2: (Color online) Coherent (solid line) and incoherent 
(dashed line) photon number in the cavity for an increasing 
number of atoms N = 100 (top), N = 1000 (middle) and N = 
10000 (bottom). Thin dotted lines show the corresponding 
quantities in the thermodynamic limit (N — > oo, calculated by 
multimode Bogoliubov theory). Parameters: uia = 1, Ac = 
—2, k = 2, Uo = and n max = 2. 



It is convenient to plot separately the coherent and in- 
coherent parts, since the N c factor gives rise to orders 
of magnitude differences. Moreover, the qualitative de- 
pendence is significantly different. The mean field re- 
flects nicely that the non-analytic point in y = y c is 
smoothed out for finite N, however, the cross-over re- 
gion becomes narrower as N is increased. In parallel, the 
number of incoherent photons increases as a power law, 
max {(a) a)} oc N T , with the exponent r = 0.41 ± 0.02. 
This result points out the sharp difference between the 
open and the closed Dicke models, since this exponent 
significantly differs from ~, which is the exact finite-size 
exponent of the closed system [24]. The bottom panel 
of the figure shows that the HFB results for N = 10000 
are practically indistinguishable from those of the Bogoli- 
ubov theory (thin dotted lines) where the limit N — » oo 
is taken. For generating the dotted line results, we ex- 
tended the model of Ref. [22] to include higher excited 
modes up to the same cutoff n max = 2 [35] . 

The very same behavior is reflected in the populations 
of the atomic modes, shown in Fig. 3 for the mean field 



FIG. 3: (Color online) The condensate mean field components 
in the Fourier modes for increasing number of atoms N = 100 
(top), N = 1000 (middle) and N = 10000 (bottom). Thin 
dotted lines represent the results of the thermodynamic limit. 
Parameters are the same as in Fig. 2. 



and on Fig. 4 for the fluctuations. In addition to elucidat- 
ing the formation of criticality, this plot reveals quantita- 
tively the negligible effect of the higher excited modes on 
the threshold region. The population I72I 2 in the mode 
cos(2fcir) starts to increase only quadratically and only 
from the critical point. It follows then that the critical 
behavior of the coupled BEC-cavity system is equivalent 
to that of the Dicke model which corresponds to the two- 
mode approximation. 

In Fig. 5 we present the spectrum of fluctuations for 
three different atom numbers. The spectrum is made 
of complex eigenfrcqucncics with non- vanishing, negative 
imaginary parts due to the cavity loss rate k. By con- 
struction, it possesses the symmetry uij <-> —ui* , hence we 
plot only half of the spectrum with positive real parts. 
The spectra is labeled according to the y — value 
of the real part, that is, the curves starting from wr, 
—Ac = 2ljr, and Aujr correspond to the mode C\ (0J1), 
the cavity field mode a (w cav ), and the motional mode C2 
(cj 2 ), respectively. These plots reveal how close the "soft" 
mode's frequency u>\ approaches zero near the critical 
point. In the thermodynamic limit (thin dotted lines), 
its real part is zero in a finite interval around y c , that is 
a distinctive feature of the open-system phase transition 
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FIG. 4: (Color online) Incoherent excitation in the Fourier 
modes of the atomic motion above the condensate mean field 
for increasing number of atoms N = 100 (top), N = 1000 
(middle) and N = 10000 (bottom). Thin dotted lines rep- 
resent the results of a multimode Bogoliubov theory of the 
thermodynamic limit. Parameters are the same as in Fig. 2. 



[22]. Besides this, the imaginary parts reflect that the 
motional modes mix with the lossy photon mode. Follow- 
ing the symmetry wi o Uj = — cj*, at the two end points 
of the interval where the real parts of lu 1 vanish, their 
imaginary parts bifurcates (only the upper branch is plot- 
ted) . The critical point is reached when the upper branch 
hits zero. The finite-size scaling is manifested in the van- 
ishing of the imaginary part, min.{|Imwi|} oc N~ e . We 
find the exponent e = 0.44±0.02. Note that Imux defines 
the correlation time of the system. The admixture of the 
mode £2 with the field mode is very small and the cavity 
decay affects this mode only marginally, what supports 
that the truncation of the Fourier set of modes, which is 
a basic ingredient in many papers, is well justified. 

Let us analyze the characteristic magnitude of finite- 
size corrections to the thermodynamic limit results. The 
finite-size effect is due to the back coupling of the second- 
order correlations of the modes a and c.\ into the mean 
field equations. It follows from (28) that the eigenfrc- 
quencies appear in the denominator, whereas k stands in 
the numerator from the diffusion matrix. The significant 
peaks of the correlation functions near the critical point, 
exhibited in Figs. 2 and 4, are due to the soft mode fre- 



FIG. 5: (Color online) Spectrum of excitations: real part (left 
column) and imaginary part (right column) for N = 100 (up- 
per row), N = 1000 (middle row) and N = 10000 (bottom 
row). Thin dotted lines show the results of the thermody- 
namic limit. Parameters are the same as in Fig. 2. 



quency, uix, that tends to vanish in the critical point. As 
the HFB terms in Eqs. (10) and (13) are multiplied by 
the condensate atom number 1/N C ~ 1/iV, the order of 
magnitude characterizing the back coupling is 



iVIm{u;i} ' v ; 

At the vicinity of the critical point, the back coupling 
scales with exponent e — 1 < 0, thus the HFB corrections 
vanish for N — > oo. It follows then that the thermody- 
namic limit of the HFB theory renders the results of the 
Bogoliubov theory. 



B. Fano factor 

In the following we compare the HFB results to those 
of an exact numerical calculation with n max = 1. This 
latter is performed by the Monte-Carlo Wave Function 
method using the general purpose C++QED package 
[36, 37] and gives numerically exact solutions for the 
equations of motion in (5). We run the simulations for 
long enough time to reach the steady-state. 

The numerical method yields the mean photon num- 
ber, for example, but it is not separated to the contri- 
butions of the mean field and that of the fluctuations. 
Therefore we need to find a quantity, other than the 
second-order noise correlations, which can characterize 
the criticality and can be assessed both in the HFB and 
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numerically. This is the Fano factor 

((at a) 2 ) - a) 2 

F =- — m — < < 32 > 

which contains the variance of the field amplitude, how- 
ever, the quantum noise part is not suppressed by the 
mean field, this latter being magnified by the factor of 
N c . For a coherent state the Fano-factor is 1 regard- 
less of the amplitude. Deviation from unity thus reveals 
the increasing role of the fluctuations with respect to the 
mean, so that the appearance of a peak at the critical 
point can be expected. The scaling of the peak height 
with the atom number can be evaluated. 

The Fano factor within the HFB approximation is 

. + . a* 2 H + a 2 (atat) / 1 \ 

This is compared to the results of the numerical calcula- 
tion in Fig. 6 as a function of the pump strength and for 
various atom numbers. In accordance with our expecta- 
tion, well-behaved peaks appear for finite N. On increas- 
ing the atom number, the peaks get narrower and higher, 
indicating how the system tends to the singular behav- 
ior of the thermodynamic limit. Note that there is some 



5 - 




0.4 0.6 0.8 1 1.2 1.4 1.6 



FIG. 6: (Color online) The Fano factor as a function of pump 
strength. Lines without line points show the HFB approxima- 
tion, lines with line points correspond to the exact numerical 
results. The peaks in increasing order correspond to atom 
numbers N = 10, 50, 100 and 200. Parameters are the same 
as in Fig. 2. 

difference in the height and the position of the corre- 
sponding peaks in the two approaches. The HFB leads to 
a peak position below, whereas the numerical simulation 
gives it above the critical point. The convergence proper- 
ties are, however, amazingly similar. In Fig. 7 the power 
law finite-size scaling of the peak height with the number 
of atoms is presented on a log-log scale. The exponent 
obtained from the fit on the HFB results (blue squares) 
0.39 ± 0.02 agrees well with the exponent acquired from 
the numerical approach (red circles) 0.40 ± 0.01. 
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FIG. 7: (Color online) Critical exponents of the finite- 
size scaling of the peak height in the Fano-factor. Solid 
line shows a linear fit on the calculated points. The 
HFB calculation is performed for atom numbers N = 
10, 100, 1000, 10000, 100000 (the fit discards the points 
N = 10 and N = 100), the numerical simulation is limited 
to N = 10, 50, 100, 200, 400 (the fit is without N = 10). 
Parameters are the same as in Fig. 2. 

C. Scaling ansatz and finite-size exponents 

To discuss the results further and to achieve a deeper 
understanding let us first review shortly the scaling re- 
lations found in the Bogoliubov theory [22], i.e., for an 
infinite system. The critical behavior of the open-system 
Dicke model follows from the vanishing of the character- 
istic frequency u>i . Interestingly the real part vanishes in 
an extended interval of y/y c and the imaginary part that 
is responsible for the stability bifurcates. The critical 
point is reached when the upper branch of the imaginary 
part reaches zero. In the infinite system Im uji ~ \y c ~ y\ i 
with the critical exponent £ = 1. 

From Eq. (28) one can deduce that the critical part of 
the second order correlation functions are dominated by 
the soft mode, that is {v^v v ) ~ (Im Wi) _1 ~ \y c — 
This critical behavior is in contrast to the case of the 
closed system Dicke model, where the soft mode only has 
a real part that vanishes as \y c — y\ x ^ 2 and consequently 
the second order correlations scale as \y c — j/|~ 1//2 . 

To demonstrate the scaling properties to the finite size 
system we assume a scaling ansatz for the imaginary part 
of the soft mode frequency, namely we assume that it 
depends on a combination of the two variables y = (y c — 
y)/y c and N rather than independently of the two. The 
most general form can be written as: 

Im Wa ~ <p(N'y), (34) 

with £ and e are yet undetermined exponents and ip(x) 
is a dimensionless function whose asymptotic values are 
fixed by the following argument. In the infinite system 
N — > oo and with y fixed the argument x — > oo and Im u\ 
has to give back the scaling in the thermodynamic limit, 
i.e., £ = 1. Therefore ]im a ._ >00 <f(x) = const. Similarly 
for the finite system and exactly y = the frequency has 
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Quantity 


Prop. 


Exp. 




max. 


0.41 


1 2/max | 


-0.39 




max. 


0.39 


1 2/max ?/c | 


-0.41 


IM"i}| 


min. 


-0.44 


1 2/min £/c | 


-0.39 


F 


max. 


0.39 


| |/max 


-0.39 



TABLE I: Summary of the finite-size exponents obtained from 
the HFB theory. A physical quantity $ scale with exponent 
r according to max{$} oc N T . The estimated errors of the 
exponents are below 0.02. 



to stay nonzero, therefore linx^o f[x) = x~ 1 . Conse- 
quently the minimum value scales as Im u)\ ~ N~ e . For 
finite N the location of the minimum value of Im u)\ is 
also shifted from zero to y*. Its scaling can be extracted 
by assuming that tp(x) has a minimum value at a finite 
x*, therefore y* — x*N~ e . 

The second order correlators (a^a), (cjci), and also F 
can be separated to a nonuniversal finite and to a critical 
part. The critical part, denoted by $, is proportional to 
(Im wi) -1 , as seen from Eq. (28). Therefore all of these 
quantities scale according to 



T 1 <p(N e y) 



(35) 



It follows then that for an infinite system all these cor- 
relation functions diverge around the critical point with 
the same exponent £ = 1. Moreover for a finite system 
all has the same finite size scaling exponent e which can 
be read from the scaling of the heights and also by the 
location of the corresponding peaks. 

The scaling exponents of the HFB theory are collected 
in Table I. It can be seen that e « 0.4 in good agreement 
for all the correlation functions both from the heights 
and also from the locations of the peaks. 



D. Atom- field entanglement 

The cavity field mode and the atomic motional modes 
become entangled in the steady state of the system. In 
the thermodynamic limit, the Bogoliubov approach leads 
to a non-vanishing finite logarithmic negativity even if 
the system is exposed to dissipation [22]. Note that the 
ground state of the closed Dicke model is the two-mode 
squeezed state of which the logarithmic negativity di- 
verges in the critical point [23]. The damping in the 
open Dicke model regularizes then the divergence, how- 
ever, the logarithmic negativity is still a non-analytic 
function of y at the critical value y c . In the following 
we present the finite-size corrections to the steady-state 
entanglement. 



The logarithmic negativity is an entanglement measure 
defined as [38] 



£Mp) = ln(||p T *||) 



(36) 



where p A means the partial transpose of the density 
operator with respect to the subsystem A and the trace 

norm is used, i.e., ||0|| =Tr{VOWy The bipartition, 

in our case, is such that subsystem A is the photon mode 
and all the atomic motional modes form the complemen- 
tary subsystem. In the HFB theory, the steady state is a 
Gaussian state for which the logarithmic negativity can 
be calculated from the symmetrically ordered correlation 
matrix [39, 40], 



C i,h = \ ({& . Qk}) , 



(37) 



where {. , .} denotes the anticommutator, and qj are 
the quadratures of the fluctuation operators, qQ — x — 
a), ?i = P = ^75 (a 1 — a), q^a = X 

■ h) , and o 2 ,-4.i = Pi = 4? ( bl - h ) . The el 



i(a'+a), qi = p = ^(a'-a), q 2j 

75 $ + h) . ^d g 2i + x = P S = -f 2 (fet - ~bj] 
cments of the correlation matrix C are readily obtained 
from the second order correlation functions of Eq. (28) 
that have been derived in the course of the iteration al- 
gorithm. It can be easily shown that the correlation ma- 
trix of the partially transposed state C Ta is obtained by 
means of the substitution p — > —p in the corresponding 
correlation matrix elements. The logarithmic negativity 
is then expressed by the symplectic eigenvalues i>j of C Ta 
[391 as 




v i < 2 



(38a) 
(38b) 



We plot the logarithmic negativity as a function of 
the pump strength for three different atom numbers 
N = 100,1000 and 10000 in Fig. 8. The atom-cavity en- 
tanglement is characterized by a broad peak with a max- 
imum near the critical point, however, the peak height 
does not scale with the number of atoms N. It reaches 
the maximum of the thermodynamic limit for all finite 
TV shown. On increasing N, the shape of the logarithmic 
negativity function tends to reproduce the singular be- 
havior with a jump of the derivative at y = y c , found in 
the thermodynamic limit. 



V. CONCLUSION 

Ultracold atoms loaded into the small volume of a high- 
finesse resonator allow for studying quantum many-body 
physics and correlated state of matter originating from 
a peculiar, long-range interaction between the particles. 
In this paper we made an important step forward in the 
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FIG. 8: The logarithmic negativity as a function of the pump 
strength. The atom numbers corresponding to the curves are 
N = 100 (top), N = 1000 (middle) and N = 10000 (bottom). 
Parameters: ujr = 1, Ac — —2, k — 2, Uo — s n max = 2. 



exploration of the self-organization process of atoms in a 
cavity, which is a quantum phase transition manifesting 
clearly the collective behavior of a quantum gas. Go- 
ing beyond the leading-order Bogoliubov-type descrip- 
tion, here we developed the Hartree-Fock theory to treat 
the long-range interaction, and in particular, the steady- 
state of the damped-driven open system. As a main merit 
of the theory, it renders the finite-size corrections to the 
bare mean-field theory, regularizes the singularity in the 
critical point. Finally, we were able to derive the finite- 
size scaling of the of the divergent correlation functions. 
The power law has been found together with the criti- 
cal exponent which is in a very good agreement with the 
outcome of an exact numerical simulation. 
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Appendix A: Heisenberg-equations after the mean-field decomposition 



For completeness, the general equations of motion (5) are written out explicitly here by using the mean-field 
decomposition (7). The terms are grouped according to the different powers of the fluctuation operators involved. 
The equation of motion for the cavity field mode is 



f-A c + M7 f M (2) 7 - i, 

(-A c + 7i7 t M (2) 7 - in) a + ua (Vm (2) c + c f M (2) 7) + -y (Vm (1) c + c f M (1) 7) + i|(t) 
cWc a + u (&MW~t + 7 t M (2) c) ~a + ]- y c^M {1 



0c 



1 



u c f M (2) c 5 



(Al) 
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whereas for the matter wave field modes 



N c i—j + i—c=^N c 
at at 



lor M (0) + -y (a* + a) M (1) + mq*«M (21 - fjl ) 7 



lorM^ + ]-y (a* + a) M (1) + ua*a M (2) - fil ) c + ( ]-y (a 1 + a) M (1) + u (a* a + a) a) M (2) ) 7 



fNr. 



~-y (& f + a) M (1) + u (a* a + a^ a) M (2) ^ c + m* aM (2) 7 



1 



ataM (2) c 



(A2) 



The mean field equations are obtained exactly by taking 
the quantum average. Then, by subtracting the mean 
field equations form the above ones, the equations of mo- 
tion for the fluctuations a and c can be deduced. 

The Hartree-Fock-Bogoliubov method relies on the fol- 
lowing approximations, in (Al), 



(gt M (2) 7 
and in (A2), 



fat 



a) M (1) c 
5 f a] M (2) c r 



5 f M (2) c w <5 f 5) M (2) c + S f M (2) (5 c) + M (2) (5 f c) 



c t 5^)M (2) 7 + 7 t M (2) (ca> 
c f a\ M (2) c + c f M (2) (ca) 



h (Sc)) 
a M (2) (a f c) 

r( 2 ) /st?\« 



m« ((ate) 

M (2) (5c) H 
a'a « (a^ 5) 

(2)?a.stA/r( 2 ) In 7\ 



These substitutions lead to the model presented in 
Sec. III. 



To be specific, the first and third rows of the right- 
hand sides give the mean field equation (10), and the 
second and fourth rows give the fluctuation dynamics in 
Eq. (13). 



[1] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 

80, 885 (2008). 
[2] S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. 

Phys. 80, 1215 (2008). 
[3] D. Jaksch and P. Zoller, Annals of Physics 315, 52 

(2005). 

[4] M. Greiner, O. Mandel, T. Esslinger, T. W. Hansen, and 

I. Bloch, Nature 415, 39 (2002). 
[5] Y. Colombe, T. Steinmetz, G. Dubois, F. Linke, 

D. Hunger, and J. Reichel, Nature 450, 272 (2007). 
[6] F. Brennecke, T. Donner, S. Ritter, T. Bourdel, M. Kohl, 

and T. Esslinger, Nature 450, 268 (2007). 
[7] S. Slama, S. Bux, G. Krenz, C. Zimmermann, and P. W. 

Courteille, Phys. Rev. Lett. 98, 053603 (2007). 
[8] K. W. Murch, K. L. Moore, S. Gupta, and D. M. 

Stamper-Kurn, Nature Physics 4, 561 (2008). 
[9] M. Wolke, J. Klinner, H. Kessler, and A. Hemmerich, 

Science (2012). 

[10] P. Miinstermann, T. Fischer, P. Maunz, P. W. H. Pinkse, 



and G. Rempe, Phys. Rev. Lett. 84, 4068 (2000). 
[11] J. K. Asboth, P. Domokos, and H. Ritsch, Phys. Rev. A 

70, 013414 (2004). 
[12] R. H. Dicke, Phys. Rev. 93, 99 (1954). 
[13] K. Hepp and E. H. Lieb, Phys. Rev. A 8, 2517 (1973). 
[14] M. Hillery and L. D. Mlodinow, Phys. Rev. A 31, 797 

(1985). 

[15] C. Emary and T. Brandes, Phys. Rev. E 67, 066203 
(2003). 

[16] K. Baumann, C. Guerlin, F. Brennecke, and T. Esslinger, 
Nature 464, 1301 (2010). 

[17] D. Nagy, G. Konya, G. Szirmai, and P. Domokos, Phys. 
Rev. Lett. 104, 130401 (2010). 

[18] K. Baumann, R. Mottl, F. Brennecke, and T. Esslinger, 
Physical Review Letters 107, 140402+ (2011). 

[19] R. Mottl, F. Brennecke, K. Baumann, R. Landig, T. Don- 
ner, and T. Esslinger, Science (2012). 

[20] P. Domokos and H. Ritsch, Phys. Rev. Lett. 89, 253003 
(2002). 



13 



[21] A. T. Black, H. W. Chan, and V. Vuletic, Phys. Rev. 

Lett. 91, 203001 (2003). 
[22] D. Nagy, G. Szirmai, and P. Domokos, Physical Review 

A 84, 043637+ (2011). 
[23] N. Lambert, C. Emary, and T. Brandes, Phys. Rev. Lett. 

92, 073602 (2004). 
[24] J. Vidal and S. Dusuel, Europhys. Lett. p. 817 (2006). 
[25] S. Dusuel and J. Vidal, Phys. Rev. B 71, 224420 (2005). 
[26] G. Liberti, F. Piperno, and F. Plastina, Phys. Rev. A 81, 

013818 (2010). 

[27] J. Reslen, L. Quiroga, and N. F. Johnson, Europhys. 

Lett. 69, 8 (2005). 
[28] R. Botet and R. Jullien, Phys. Rev. B 28, 3955 (1983). 
[29] D. Nagy, G. Szirmai, and P. Domokos, Eur. Phys. J. D 

48, 127 (2008). 

[30] J. Keeling, M. J. Bhaseen, and B. D. Simons, Phys. Rev. 
Lett. 105, 043001 (2010). 

[31] S. Gopalakrishnan, B. L. Lev, and P. M. Goldbart, Na- 
ture Physics 5, 845 (2009). 

[32] S. F. Vidal, G. De Chiara, J. Larson, and G. Morigi, 
Phys. Rev. A 81, 043407 (2010). 

[33] C. Cohen- Tannoudji, J. Dupont-Roc, and G. Grynberg, 
Atom-Photon Interactions. Basic Processes and Applica- 
tions (Wiley- VCH, Weinheim, 2004). 

[34] A. Griffin, Physical Review B 53, 9341 (1996). 

[35] G. Konya, G. Szirmai, and P. Domokos, The European 
Physical Journal D - Atomic, Molecular, Optical and 
Plasma Physics 65, 33 (2011). 



[36] A. Vukics and H. Ritsch, The European Physical Journal 
D - Atomic, Molecular, Optical and Plasma Physics 44, 
585 (2007). 

[37] A. Vukics, Computer Physics Communications 183, 1381 
(2012). 

[38] M. B. Plenio and S. Virmani, Quantum Info. Comput. 7, 
1 (2007). 

[39] G. Vidal and R. F. Werner, Phys. Rev. A 65, 032314 
(2002). 

[40] G. Adesso, A. Serafini, and F. Illuminati, Phys. Rev. A 
70, 022318 (2004). 

[41] This assumption will only be required later, when we 
consider the steady state of the system, but it simplifies 
the calculation without losing interesting effects. 

[42] Since the condensate breaks the U(l) symmetry of the 
model, the Goldstone-theorem states that the spectrum 
of the system must contain a zero mode. This zero mode 
is bo and it describes the particle number and phase fluc- 
tuations of the condensate. But in the HFB approxima- 
tion, the frequency of the bo mode is not strictly zero: 
this is a well-known weakness of the HFB approxima- 
tion. This frequency is very small, but couples also to 
a and a\ which renders the problem to be numerically 
unstable. To avoid it, we set to zero the frequency and 
the couplings related to the mode bo in the matrix F by 
hand. The algorithm proved to be numerically stable. 



